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Radar  imaging  and  parameter  estimation  can  be  modeled  as  an  inverse  problem:  y  =  R(f)  +  w,  i.e.,  from  noisy  measurement 
vector  y  to  reconstruct  or  estimate  the  target  parameter  vector  f,  which  could  be  the  spatial  reflectivity  coefficient  vector  or  other 
target  properties.  The  operator  R  (linear  or  nonlinear)  depends  on  many  factors  such  as  the  waveform  transmission  strategy, 
the  data  collection  geometry,  the  medium  property,  and  the  underlying  wave  propagation  models.  A  fundamental  question  that 
this  project  aims  to  address  is:  how  to  explore  novel  waveform  transmission,  modulation,  and  statistical  computing  techniques 
to  achieve  superior  imaging  and  parameter  estimation. 

The  main  accomplishments  of  the  project  during  this  reporting  period  are  (1)  the  development  of  multi-parameter  estimation 
methods  for  cognitive  radar  under  compound  Gaussian  clutter  using  variational  Bayesian  inference;  2)  the  development  of 
range  and  Doppler  estimation  method  using  weighted  OFDM  modulation,  and  3)  the  development  of  novel  microwave  nonlinear 
image  reconstruction  algorithm  using  waveform  encoding  schemes.  A  summary  of  the  accomplishments  is  provided  below. 
Detailed  technical  descriptions  are  provided  in  the  attached  publications  and  manuscripts. 

Accomplishment  1 :  Development  of  Multi-Parameter  Estimation  Method  for  Cognitive  Radar  Under  Compound  Gaussian  Clutter 
Using  Variational  Bayesian  Inference 

Cognitive  radar  has  been  proposed  as  a  fully  adaptive  radar  transmission  and  reception  system  [1].  In  cognitive  radar,  both  the 
transmitter  and  the  receiver  parameters  are  estimated  and  updated  by  learning  from  the  unknown  environment,  forming  a  belief 
on  what  is  learned,  and  propagating  this  belief  by  Bayesian  inference.  From  the  parameter  estimation  perspective,  the  Bayesian 
approach  enables  inclusion  of  prior  information  (knowledge)  of  radar  target  and  clutter  by  estimating  the  posterior  density  of  the 
unknown  parameters.  The  estimation  is  optimal  in  the  sense  of  minimizing  the  Bayesian  mean  squared  error  (MSE).  Typically, 
the  full  joint  probability  density  function  (pdf)  of  all  the  parameters  of  interest  including  the  nuisance  parameters  is  considered. 
However,  in  the  case  of  high-dimensional  multi-variate  integration  of  Bayesian  posterior  density,  the  calculation  of  the  posterior 
pdf  and  its  marginal  can  be  computationally  prohibitive  and  tractable  analytical  solutions  are  often  not  available.  Furthermore, 
the  estimation  accuracy  is  directly  related  to  the  number  of  data  samples  in  the  Bayesian  estimator.  However,  in  many  radar 
applications,  the  number  of  available  data  samples  is  limited.  These  computational  challenges  and  limitations  must  be 
addressed  to  develop  next  generation  cognitive  radar. 

In  this  work,  we  developed  a  variational  Bayesian  (VB)  based  method  for  parameter  estimation.  Variational  Bayesian  aims  to 
minimize  free  energy  [2,  3],  which  is  equivalent  to  minimizing  the  Kullback-Liebler  divergence  between  the  true  posterior  density 
and  an  approximation  density  of  the  parameters  to  be  estimated.  As  a  result  of  this  functional  optimization  for  density 
estimation,  the  marginal  VB  posterior  density  has  an  explicit  functional  structure,  thus  leading  to  closed  form  solutions  [4].  In  this 
work,  we  focus  on  multiple  parameter  estimation  of  target  under  compound  Gaussian  clutter  in  the  context  of  cognitive  radar  as 
an  extension  of  our  prior  work  on  single  parameter  estimation  [2].  The  compound  Gaussian  clutter  model  is  used  in  high- 
resolution  and  low-grazing-angle  radar  to  characterize  random  and  non-stationary  sea  clutter.  The  Bayesian  estimator  must 
consider  a  multi-parameter  estimation  problem  by  which  the  parameters  in  the  compound-Gaussian  model,  the  radar  target 
response,  as  well  as  other  nuisance  parameters  are  estimated.  We  compare  the  performance  of  the  proposed  VB  method  with 
the  expectation-maximization  (EM)  algorithm.  In  the  EM  method,  expectations  of  sufficient  statistics  are  computed  with  respect 
to  the  posterior  density  of  hidden  variables  and  then  used  to  iteratively  estimate  the  unknown  parameters  by  the  maximum 
likelihood  principle.  We  show  that  our  proposed  variational  algorithms  outperform  the  EM  method  particularly  when  estimating 
parameters  that  follow  non-Gaussian  nonlinear  models  in  Bayesian  inference.  Hence,  the  proposed  variational  algorithms 
provide  appealing  computational  advantages  for  cognitive  radar.  The  results  have  been  presented  at  the  IEEE  Conference  on 
Acoustics,  Speech,  and  Signal  Processing  in  April  2015  (see  Ref.  [5]).  A  more  detailed  discussion  of  the  VB  algorithm  has  been 
included  in  a  journal  version  of  the  paper  (see  Ref.  [6]). 


Accomplishment  2:  Development  of  range  and  Doppler  estimation  using  weighted  OFDM  modulation  for  radar  targets. 

Estimation  of  range  and  velocity  of  targets  is  an  important  problem  in  radar  applications.  Orthogonal  frequency  division 
multiplexing  (OFDM)  offers  robust  system  performance  in  rich  multipath  fading  environment  [7].  Challenges  with  range  and 
velocity  estimation  of  radar  targets  are  the  need  for  high  range  and  Doppler  resolutions  and  resolving  phase  wrapping.  A  high 
range  resolution  requires  a  high  signal  bandwidth  while  a  high  Doppler  resolution  requires  long  pulse  repetition  intervals.  For 
OFDM  signals,  the  range  resolution  depends  on  the  sub-carrier  spacing.  The  resolution  improves  as  the  spacing  increases. 
However,  since  the  sub-carrier  spacing  and  the  pulse  repetition  interval  are  inversely  related,  there  exist  trade-offs  for  achieving 
high  resolutions  in  both  range  and  Doppler.  Next,  for  fast  moving  targets  such  as  missiles  and  airplanes,  phase  wrapping  may 
appear  due  to  the  2  pi  modulo  folding,  i.e.,  the  actual  Doppler  frequency  is  likely  to  be  greater  than  the  corresponding  sampling 
frequency.  As  a  result,  discrete  Fourier  transform  (DFT)  based  sinusoidal  estimation  methods  can  only  extract  the  remainder 
after  2  pi  modulo  folding.  The  integer  multiple  of  the  sampling  frequency  that  is  lost  due  to  phase  wrapping  is  known  as  the 
folding  error.  The  time  delay  estimation  corresponding  to  the  range  is  also  subject  to  a  similar  phase  wrapping  problem. 
Furthermore,  phase  based  parameter  estimation  methods  are  generally  sensitive  to  phase  noise.  Hence,  for  the  OFDM 
signaling  scheme,  estimation  method  that  resolves  both  the  remainder  and  the  folding  error  for  range  and  velocity  is  needed. 


In  this  work,  we  extend  our  prior  work  on  radar  OFDM  modulation  [8],  [9]  and  addresses  the  problem  of  estimation  of  range  and 
velocity  of  a  radar  target  and  OFDM  waveform  design.  The  contribution  of  our  work  is  threefold.  First,  we  propose  a  two  stage 
estimation  method  to  estimate  the  fractional  components  and  the  folding  integers  of  range  and  velocity,  respectively,  using 
OFDM  waveforms.  In  the  first  stage,  we  employ  the  maximum  likelihood  approach  to  estimating  the  fractional  components  of 
the  target  parameters.  In  the  second  stage,  the  robust  Chinese  Remainder  theorem  is  utilized  to  extract  folding  integers  in 
parameter  estimation.  For  this  purpose  we  present  a  design  of  variable  frequency  step  in  the  OFDM  signaling  for  facilitating 
robust  phase  unwrapping.  Second,  we  derive  the  Cramer-Rao  lower  bounds  for  the  maximum  likelihood  estimator  for  range  and 
Doppler  and  the  total  error  variance  for  the  developed  two-stage  estimator.  Third,  the  weights  of  the  OFDM  symbols  are 
designed  by  optimizing  the  error  bound  of  the  estimator  subject  to  the  constraints  on  the  peak  to  average  power  ratio  (PARP) 
and  the  total  transmission  energy.  Numerical  simulations  show  that  the  weighted  OFDM  scheme  improves  the  Cramer-Rao 
bound  on  the  range  estimation  accuracy  while  controlling  the  maximum  level  of  PAPR.  Hence,  the  weighted  OFDM  modulation 
method  provides  a  flexible  modulation  mechanism  for  radar  with  an  improved  range  estimation  accuracy.  A  provisional  patent 
on  this  technology  has  been  filed  to  the  US  Patent  and  Trademark  Office  [10]. 


Accomplishment  3:  Waveform  Encoding  for  Nonlinear  Electromagnetic  Tomographic  Imaging 

Electromagnetic  (EM)  tomographic  imaging  is  an  inverse  scattering  problem  which  has  a  wide  range  of  applications  in  medical 
imaging,  geophysical  exploration,  nondestructive  testing,  and  target  identifications.  In  EM  tomographic  imaging,  source 
antennas  transmit  EM  signals  into  a  medium  under  test  and  receive  scattering  signals.  Based  on  the  underlying  Maxwell’s 
equations,  inversion  methods  are  employed  to  reconstruct  a  spatial  distribution  of  material  parameters  such  as  the  dielectric 
permittivity  and/or  magnetic  permeability  of  the  target  and  the  surrounding  medium,  thus  turning  recorded  scattering  data  into 
images. 

The  imaging  problem  is  formulated  mathematically  as  a  nonlinear  optimization  problem  that  seeks  to  minimize  a  misfit  function 
between  the  measured  sensor  data  and  a  pre-determined  model  conditioned  on  parameters  which  are  to  be  reconstructed.  An 
iterative  gradient-descent  algorithm  is  developed  to  solve  the  inverse  operator  (also  called  the  adjoint  operator)  problem.  What  it 
means  that  the  image  (i.e.,  the  spatial  distribution  of  the  material  property  values)  to  be  reconstructed  is  to  be  updated 
iteratively  until  a  stopping  criteria  is  reached.  We  call  this  method  the  propagation  and  backpropagation  (PBP)  method.  To  be 
more  specific,  the  propagation  step  is  to  calculate  a  predicted  value  of  the  wave  field  data  from  a  forward  model,  while  the 
backpropagation  step  is  to  update  the  image  value  through  the  use  of  adjoint  method.  For  iterative  algorithms,  slow 
convergence  and  high  computational  complexity  are  the  limiting  factors  for  real-time  applications. 

In  this  work,  we  will  address  this  issue  by  means  of  dimensionality  reduction  through  multiple  source  waveform  encoding. 

Typical  full  wave  tomographic  imaging  operates  in  a  single-input  multiple-output  (SIMO)  configuration.  An  image  is 
reconstructed  from  measured  data  in  response  to  a  single  excitation  antenna  source.  The  reconstruction  process  continues  till 
all  the  sources  are  excited.  For  large-scale  imaging  such  as  seismic  exploration,  the  number  of  EM  sources  is  very  large.  Not 
only  the  computational  cost  is  high,  the  operational  expenditures  of  each  data  collection  process  is  also  significant.  In  this  case, 
multiple  source  excitation  becomes  appealing.  For  example,  in  seismic  imaging,  multiple  sources  are  excited  simultaneously  to 
form  a  so-called  supershot  to  probe  the  imaging  field  [11],  which  means  the  imaging  configuration  becomes  multiple-input 
multiple-output  (MIMO).  The  recorded  data  are  processed  to  form  an  image.  The  image  is  updated  when  new  measurement 
data  is  available.  This  procedure  is  repeated  until  the  image  converges  or  a  predetermined  stopping  criterion  is  met.  However, 
multiple  wave  simultaneous  excitation  induces  cross-talk  noise  due  to  wave  interference,  which,  if  not  treated,  will  cause  image 
artifacts.  Therefore,  signal  processing  techniques  such  as  waveform  encoding  are  needed  to  mitigate  cross-talk  noise  in  the 
image  in  order  to  achieve  high  quality  imaging  while  reducing  the  computational  complexity.  The  contribution  of  our  work  is 
threefold.  First,  we  develop  three  different  waveform  encoding  techniques,  i.e.,  random  phase  encoding,  waveform  delay 
encoding,  and  uniform  weight  encoding  for  the  full-wave  EM  imaging  problem.  We  show  that  the  random  phase  encoding 
method  results  in  constant-envelop  waveforms  and  produces  the  best  performance  in  terms  of  convergence.  Second,  this  paper 
extends  our  early  work  on  microwave  imaging  in  a  SIMO  configuration  [12].  We  show  that  using  simultaneous  sources  made  of 
superposition  of  encoded  sources  is  able  to  accelerate  iterative  algorithms  for  electromagnetic  full-wave  inversion,  thus 
demonstrating  the  effectiveness  of  waveform  encoding,  a  common  signal  processing  technique,  to  improve  computational 
efficiency  of  classic  nonlinear  inverse  problems.  Third,  although  waveform  encoding  techniques  have  been  studied  for 
acoustical  wave  imaging  (see  our  early  work  on  MIMO  ultrasonic  imaging  [13],  [14]),  there  is  still  a  lack  of  research  for  EM  full- 
wave  imaging  in  applications  where  the  use  of  EM  wave  energy  is  critical.  In  this  work,  we  develop  iterative  algorithms  that 
solve  time  domain  Maxwell’s  equations  with  coupled  electric  fields  and  magnetic  fields  using  waveform  encoded  excitations  and 
demonstrate  faster  convergence  compared  with  our  prior  SIMO  imaging  algorithm  in  [12].  The  results  have  been  submitted  to 
the  2015  IEEE  Global  Conference  on  Signal  &  Information  Processing  (see  Ref.  [15]) 


In  summary,  during  this  project  period,  we  have  addressed  successfully  the  overarching  technical  problem  of  the  project  in 
three  aspects,  i.e.,  multi-parameter  estimation  in  cognitive  radar,  OFDM  modulation  design  for  range  and  Doppler  estimation, 
and  waveform  encoding  for  microwave  imaging.  Our  research  has  led  to  two  top-tier  IEEE  journal  publications  (see  Ref.  [16], 


[17]),  one  peer-reviewed  primary  conference  publication  (see  Ref.  [5]),  one  conference  submission  (see  Ref.  [15]),  and  three 
other  journal  submissions  (see  Ref.  [6],  [18],  [19]).  A  provisional  patent  is  also  filed  to  the  US  Patent  and  Trademark  Office  in 
May  2015  (see  Ref.  [10]). 
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Abstract 

Radar  imaging  and  parameter  estimation  can  be  modeled  as  an  inverse  problem:  y  =  R  •  f  +  w,  i.e., 
from  noisy  measurement  vector  y  to  reconstruct  or  estimate  the  target  parameter  vector  f ,  which 
could  be  the  spatial  reflectivity  coefficient  vector  or  other  target  properties.  The  operator  R  (linear 
or  nonlinear)  depends  on  many  factors  such  as  the  waveform  transmission  strategy,  the  data  collec¬ 
tion  geometry,  the  medium  property,  and  the  underlying  wave  propagation  models.  A  fundamental 
question  that  this  project  aims  to  address  is:  how  to  explore  novel  waveform  transmission,  modula¬ 
tion,  a  priori  knowledge  of  target,  and  statistical  computing  techniques  to  achieve  superior  imaging 
and  parameter  estimation. 

The  main  accomplishments  of  the  project  during  this  reporting  period  arc:  (1)  The  develop¬ 
ment  of  multi-parameter  estimation  methods  for  cognitive  radar  under  compound  Gaussian  clutter 
using  variational  Bayesian  inference;  2)  The  development  of  range  and  Doppler  estimation  method 
using  weighted  OFDM  modulation;  and  3)  The  development  of  novel  microwave  nonlinear  image 
reconstruction  algorithm  using  waveform  encoding  schemes. 
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Chapter  1 

Multi-Parameter  Estimation  by 
Variational  Bayesian 

1.1  Introduction 

Cognitive  radar  has  been  proposed  as  a  fully  adaptive  radar  transmission  and  reception  system 
in  [Hay kin,  2006].  In  cognitive  radar,  both  the  transmitter  and  the  receiver  parameters  arc  esti¬ 
mated  and  updated  by  learning  from  the  unknown  environment,  forming  a  belief  on  what  is  learned, 
and  propagating  this  belief  by  Bayesian  inference.  From  the  parameter  estimation  perspective,  the 
Bayesian  approach  enables  inclusion  of  prior  information  (knowledge)  of  radar  target  and  clutter 
by  estimating  the  posterior  density  of  the  unknown  parameters.  The  estimation  is  optimal  in  the 
sense  of  minimizing  the  Bayesian  mean  squared  error  (MSE).  Typically,  the  full  joint  probability 
density  function  (pdf)  of  all  the  parameters  of  interest  including  the  nuisance  parameters  is  con¬ 
sidered.  However,  in  the  case  of  high-dimensional  multi-variate  integration  of  Bayesian  posterior 
density,  the  calculation  of  the  posterior  pdf  and  its  marginal  can  be  computationally  prohibitive 
and  tractable  analytical  solutions  arc  often  not  available.  Furthermore,  the  estimation  accuracy  is 
directly  related  to  the  number  of  data  samples  in  the  Bayesian  estimator.  However,  in  many  radar 
applications,  the  number  of  available  data  samples  is  limited.  These  computational  challenges  and 
limitations  must  be  addressed  to  develop  cognitive  radar. 

In  [Turlapaty  and  Jin,  2013],  we  proposed  a  variational  Bayesian  (VB)  based  method  for  param¬ 
eter  estimation  and  waveform  design  where  a  single  parameter  estimation  problem  is  considered. 
Variational  Bayesian  aims  to  minimize  free  energy  (FE)  [Friston,  2010,  Turlapaty  and  Jin,  2013], 
which  is  equivalent  to  minimizing  the  Kullback-Fiebler  divergence  between  the  true  density  and 
an  approximation  density.  As  a  result  of  this  functional  optimization  for  density  estimation,  the 
marginal  VB  posterior  density  has  an  explicit  functional  structure,  thus  leading  to  closed  form  so¬ 
lutions  [Smidl  and  Quinn,  2008].  This  work  extends  our  prior  work  in  [Turlapaty  and  Jin,  2013]  to 
multiple  parameter  estimation  in  the  context  of  cognitive  radar.  In  adaptive  radar  detection,  estimat¬ 
ing  the  clutter  covariance  matrix  is  a  very  important  task  since  the  detection  performance  depends  di¬ 
rectly  on  the  accuracy  of  the  estimate.  For  example,  in  high-resolution  and  low-grazing-angle  radar, 
only  a  small  sea  surface  area  is  illuminated  by  a  narrow  radar  beam.  The  sea  clutter  due  to  reflection 
from  the  small  patch  of  sea  surface  is  random  and  non- stationary  [Sangston  and  Gerlach,  1994], 
which  is  commonly  modeled  as  a  compound-Gaussian  distribution  to  characterize  its  heavy-tailed 
clutter  distributions  [Wright,  1968,  Greco  et  ah,  2004].  Hence,  the  Bayesian  estimator  must  con¬ 
sider  a  multi-parameter  estimation  problem  by  which  the  parameters  in  the  compound-Gaussian 
model,  the  radar  target  response,  as  well  as  other  nuisance  parameters  arc  estimated.  In  this  work. 
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we  compare  the  performance  of  the  proposed  VB  method  with  the  expectation-maximization  (EM) 
algorithm.  In  the  EM  method,  expectations  of  sufficient  statistics  are  computed  with  respect  to  the 
posterior  density  of  hidden  variables  and  then  used  to  iteratively  estimate  the  unknown  parame¬ 
ters  by  the  maximum  likelihood  principle  [Candy,  2009,  Wang  et  ah,  2006].  We  show  in  this  work 
that  the  variational  algorithms  outperform  the  EM  method  particularly  when  estimating  parameters 
that  follow  non-Gaussian  nonlinear  models  in  Bayesian  inference.  Hence,  the  proposed  variational 
algorithms  provide  appealing  computational  advantages  for  cognitive  radar. 

1.2  Problem  Formulation 

1.2.1  Radar  signal  model 

The  compound  clutter  model  is  a  product  of  two  random  processes  [Wang  et  ah,  2006,  Wright,  1968, 
Greco  et  ah,  2004], 

ipt  =  y/utwt  (11) 

where  the  speckle  wt  characterizing  the  local  scattering  and  is  modeled  as  a  zero  mean  complex 
Gaussian  (ZMCG)  process  wt  ~  CJ\f{0,  cr2).  The  component  iq  is  a  slow  changing  process  termed 
texture  that  follows  an  inverse  Gamma  distribution  ut  ~  T  _1(a),  where  the  pdf  of  iq  is 

p{ut\ a)  =  exp  \  (1.2) 

T(a)  V  utJ 

The  model  (1.1)  is  referred  to  as  K-clutter  [Gini  et  ah,  2000,  Sangston  and  Gerlach,  1994]  and  the 
parameter  a  is  known  as  the  Nakagami  parameter  [Haykin  et  ah,  2002].  Next,  we  assume  the  radar 
transmits  a  waveform  (I>/  and  the  electromagnetic  (EM)  energy  hits  a  target  with  a  complex  ampli¬ 
tude  response  x.  The  reflected  EM  energy  is  intercepted  by  the  radar  receiver.  The  radar  signal 
model  that  includes  the  clutter  from  (1.1)  is  given  by 

Ut  =  ®tx  +  V’t,  t  =  1,2,-  ■  •  ,N  (1.3) 

The  conditional  probabilistic  model  of  the  measurements  is  a  complex  Gaussian  distribution  given 
by 

yt\ut,x,a2  ~  CAT($tX,utcr2)  (1.4) 

1.2.2  Problem  of  parameter  estimation 

The  complete  hierarchical  stochastic  model,  i.e.,  the  joint  probability  density  function  of  the  mea¬ 
surements,  hidden  variables,  and  the  unknown  parameters  at  time  t  is  given  by 

p(yt,ut,x,cr2,a;$t)  =  p(yt\ut,x,a2-,<S>t)p0(x,cr2)p(ut\a)p0(a)  (1.5) 

where  p(yt\ut,  x,  a 2;  4>/ )  is  the  conditional  density  given  in  (1.4).  The  probabilistic  model  of  tex¬ 
ture  p(ut\a)  is  given  by  (1.2).  p0(x,  a2)  and  /;0(aj  are  prior  densities  of  (x.  a2)  and  a,  respectively. 
Initially,  the  unknown  parameter  vector  is  given  by  [x,  cr2,  a].  When  we  use  a  variational  estimation 
method,  the  covariance  estimate  depends  on  the  current  estimate  of  a.  However,  initially  this  mutual 
dependence  of  estimates  introduces  a  multiplicative  error  in  the  estimates  of  these  two  parameters. 
To  correct  this  error,  we  use  the  technique  of  covariance  adjustment  and  introduce  an  additional  pa¬ 
rameter  A  in  the  texture  model  (1.2)  by  redefining  the  covariance  in  terms  of  an  adjusted  covariance 
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<r2,  [Wang  et  ah,  2006],  i.e.,  the  relation  between  the  actual  covariance  and  the  adjusted  covariance 
is  a 2  =  <r2/ A.  Hence,  the  adjusted  texture  model  becomes 

p(ut\a,\)  =  (a\)a /T(a)  u^a~1e~a^Xut^  (1.6) 

The  new  augmented  parameter  vector  to  be  estimated  is 

0  =  [; x ,  a2,  a,  A]  (1.7) 


1.3  Variational  Bayesian  Estimator 

1.3.1  Background  of  variational  Bayesian  inference 

Consider  a  hierarchical  probabilistic  model 

p(  Y,  X,  6)  =  p{  Y|X,  0)p(X|0)p(0)  (1.8) 

where  Y  arc  the  measurements,  X  arc  hidden  variables,  and  0  arc  unknown  parameters.  In  the 
exact  Bayesian  approach,  the  unknown  parameters  are  determined  by  evaluating  the  joint  posterior 
density  p(X,  0|Y)  using  the  Bayes  rule 

p(X,  0|  Y)  =  p(Y ,  X,  6)/p(Y)  (1.9) 

while  the  marginal  posterior  p(6 |Y)  is  evaluated  by  integrating  X  out  from  the  joint  posterior. 
However,  in  practice  the  denominator  in  (1.9)  is  usually  theoretically  intractable  except  in  some 
special  cases.  Moreover,  in  the  case  of  multiple  parameters  in  6  even  the  numerical  integration 
is  computationally  expensive  and  time  consuming  [Smidl  and  Quinn,  2008].  To  address  this  issue, 
approximation  methods  arc  needed  to  determine  alternative  density  functions.  Variational  Bayesian 
estimation  aims  to  find  approximations  <y(X)  and  q(0)  to  the  marginal  posterior  densities  of  pa¬ 
rameters  that  minimize  the  variational  free  energy  of  the  approximate  density  and  the  joint  pdf  in 
(1.8)  [Beal  and  Ghahramani,  2003].  So  the  key  idea  is  factorization  of  c/(X.  6)  into  q(X)q(0),  thus 
separating  the  densities  of  X  and  6.  The  variational  free  energy  (FE)  for  Y,X  and  6  is  given 
by  [Turlapaty  and  Jin,  2013,  Smidl  and  Quinn,  2008,  Friston,  2010] 

F( Y,X,0)  =  -  J  q(X)q(0)  log  d*d°  (1.10) 

The  approximate  density  q(X)  of  hidden  variables  is  determined  by  variational  Bayesian  expecta¬ 
tion  (VBE) 

q(X)  =  argminF(Y,X,  6)  (1.11) 

?(x) 

where  the  optimal  solution  to  (1. 1 1)  is  given  by 

logry(X)  =  (logp(X|0))9(0)  +  (logp(Y|X,  O))q{0)  +  const.  (1.12) 

where  (•)  is  an  inner  product.  The  approximate  density  q(0)  of  the  parameters  is  determined  by 
variational  Bayesian  minimization  (VBM) 

q(6)  =  argminF(Y,  X,  0)  (1-13) 

q(0) 

where  the  optimal  solution  to  (1.13)  is  given  by  [Paisley  et  ah,  2012,  Tzikas  et  ah,  2008] 

log q(0)  =  logpo(0)  +  (logp(Y|X,  0))9(X)  +  const.  (1.14) 
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1.3.2  Variational  Bayesian  estimator 

The  new  joint  probability  density  function  depends  on  the  measurement  vector  Y  =  [yi .  •  ■  •  ,  Vn]T , 
the  vector  of  independent  hidden  variables  U  =  [u i ,  •  •  •  ,un]t,  and  the  adjusted  parameter  vector 
6,  which  takes  the  form  of 

p(yt,ut,x,a2,a,  A;T>t)  =  p(yt\ut,x,  a2;  &t)p0(x,  a2)p(ut\a,  X)p0(a,  A)  (1.15) 
The  full  data  log-likelihood  function  is  given  by 

A(Y,  U,  0)=^  log P{yt\ut,  x ,  cr2;  (ht)+^  log  p(ut\a,  A)+logp0(.x,  cr2)+logp0(a,  A)  (1.16) 

t  t 

Eqn.  (1.16)  involves  six  parameters,  which  is  very  complicated.  Next,  we  present  a  variation 
Bayesian  approach  to  approximate  it  in  the  sense  of  minimization  of  the  free-energy  (i.e.,  KL  diver¬ 
gence,  [Turlapaty  and  Jin,  2013]),  so  that  analytically  trackable  closed  form  expressions  of  pdfs  can 
be  obtained.  The  evaluation  of  (1.16)  involves  the  expectation  and  minimization  steps  (i.e..  E-step 
and  M-step)  while  using  the  VB  approach.  Due  to  space  limitations,  we  will  omit  some  details  of 
the  mathematical  derivation. 

VB  estimation  step  (VBE) 

The  first  step  in  variational  estimation  is  to  determine  the  approximate  density  of  the  hidden  vari¬ 
able  ut  given  the  measurements  yt  by  using  the  minimum  free  energy  principle.  Since  the  hidden 
variables  iit  arc  independent  across  time,  the  approximate  density  of  the  vector  U  can  be  written 
as  g(U)  =  H/V  i  q(ut)-  Us  approximate  density  will  be  obtained  by  minimizing  the  free-energy 
defined  in  (1.10),  which  can  be  re-written  as 


N 

F(Y,U,0)=-(A(Y,U,0)\  +(X>g(</K)))  ,rn 

based  upon  (1.16).  Hence,  the  optimization  problem  is  given  by 

q(ut)  =  argminT’(Y,  U,  6)  (1-17) 

<?(wt) 

Following  the  optimal  solution  in  (1.12),  we  obtain 

log q{ut)  =  (log p(ut\a,  X))q(a,\)  +  (log p(yt\ut,  x,  cr2))q(x^2)  +  const.  (1.18) 
Next,  inserting  the  pdf  (1.4)  in  (1.18),  we  have 

q(ut)  ex  u^(<a>«(«‘.A)+1)_1e-[<a/A>«(a.A)+<lw-***laj|>,(ie,„2)]/«*  (1<19) 

which  is  consistent  with  the  definition  of  an  inverse  Gamma  density  function  T_1(-),  leading  to  a 
closed  form  q(ut)  =  T_1  [ut\ cu,  djj')  with  its  parameters  being  c\j  =  (a)9(a  ^)  +  1  an<i  djj  = 

(°l/\ } q(a,\)  (I Vt  /&  ) q(x,a2)' 
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VB  minimization  step  (VBM) 

The  optimization  problem  for  estimating  the  pdf  of  the  parameter  vector  6  is  given  by 

q(x,  (T2,  a,  A)  =  arg  min  F(Y,  U,  9)  (1-20) 

q(x,a2  ,a,\) 

We  write  the  variational  posterior  density  of  the  parameter  vector  in  (1.14)  as  follows 

N 

log  <70,cr^,  a,  A)  =  logp0(a,  \)  +  logp0(x,  (7%)  +  (logp(yt\ut,  x,  crl;  <S>t))  ^ut)  + 

t=  l 
N 

Y  ( logp(«t|a,  A)>?(ut)  +  const.  (1.21) 

t=  i 

Note  that  the  terms  for  the  joint  pdf  of  (x,  a%)  are  separable  from  the  pdf  of  (a,  A),  i.e., 

q{x,  0-2,  a,  A)  =  q(x,  a2)  q(a ,  A)  (1.22) 

Next,  we  evaluate  logq(x,  aa)  and  log  q(a,  A)  from  (1.21)  and  (1.22)  to  derive  approximate  closed 
forms.  We  stai't  by  writing  the  log-function  of  pdf  of  (x,  a2)  as 

N 

\ogq(x,al)  =  logp0(x,al)  +  Y(log1/(‘Kutcrl))  (1.23) 

N 

~Y(\yt~  ,  +  COnst- 

X  '  ?(«t) 

which  leads  to  the  joint  pdf  q(x.  a2),  which  is  a  complex  Gaussian  inverse  Gamma  {CQTQ)  distri¬ 
bution 


q(x,ai)  oc  p0O,a2)(l/a2)^e-S*=il^-^-l2/^(iM>9(«t)  (1.24) 

Note  that  the  functional  form  is  nicely  preserved  if  the  prior  pdf  is  assumed  to  be  p0(x,a2)  ~ 
CQlQ(p,q,  /3,  px),  which  leads  to  an  explicit  single  CQTQ  distribution  expression  for  q(x,  a2). 
Next,  we  examine  the  logarithm  of  the  joint  pdf  of  (a,  A),  which  is  given  by  (where  C  is  a  constant) 

N 

log q(a,  A)=logp0(a,  A)+^(logp(nt|a,  A))g(ut)+  C  (1.25) 

t= l 

By  (1.6)  and  7)  =  Y.t=i(lo&ut)q(Ut),  T{  =  YZ.ii1  /  ut)  q{uty  we  obtain  the  joint  pdf  q(a,  A)  as 
follows 

fa  \Na 

q(a,  A)  (X  p0(a,  A)  2*  (T26) 

(r(a)) 

Since  this  is  not  a  known  probability  distribution,  it  can  not  provide  closed  form  solutions.  We 
can  rely  on  numerical  means  for  calculation,  or  alternatively,  we  propose  to  induce  another  layer  of 
factorization  on  the  joint  pdf  q(a,  A)  as 


q{a,  A)  =  q(a)q{ A) 


(1.27) 
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via  the  optimization  problem  to  obtain  closed  form  expressions  of  q{a)  and  q( A),  i.e., 

arg  min  F(Y,  U,  x,  a2,  a,  A)  (1-28) 

<?(«), ?Q) 

where  the  new  free-energy  quantity  is  defined  by 


F(Y,U,x,cr2,a,  A)  =  — <^A(Y,  U,  x,  a2,  a,  A)^ 


q(a)q(\)q(x,a2)q(U) 


+{q(a)q(X)q(x,  a2)q( U) 


q(a)q(\)q(x.,a2)q{U) 


As  a  result,  we  obtain  the  approximate  inverse  Gamma  distribution  for  q( A)  =  I'- 1  (A:  c\.  d\) 
and  the  approximate  Gamma  distribution  for  q(a)  =  T(ca(0),  dQ(0))  by  Lindley’s  approxima¬ 
tion  [Lindley,  1980].  Due  to  space  limitations,  we  will  omit  the  detailed  expressions  of  q(a)  and 
q( A)  in  this  work. 


1.4  Numerical  Simulations 

We  present  estimation  performance  of  the  three  estimators  by  numerical  simulations.  The  three  esti¬ 
mators  are:  (1)  VBN-PF:  Variational  Bayesian  approach  with  induced  factorization  on  q(a,  A)  and 
a  numerical  integration  for  evaluation  of  the  moments  of  a.  (2)  VBL-PF :  Variational  Bayesian  ap¬ 
proach  with  induced  factorization  on  q(a,  A)  and  the  Lindley’s  approximation  for  evaluation  of  mo¬ 
ments  of  a  [Lindley,  1980].  (3)  PXEM:  parameter  expanded  expectation  maximization  approach, 
where  q{ut)  is  updated  based  on  Bayesian  approach,  while  each  of  the  unknown  parameters  is  deter¬ 
mined  by  using  the  maximum  likelihood  principle  [Wang  et  ah,  2006].  The  key  differences  between 
the  EM  and  VB  algorithms  are  twofold.  First,  in  the  EM  method,  unknown  parameters  are  consid¬ 
ered  as  deterministic  values  and  estimated  using  the  maximum  likelihood  (ML)  method  whereas 
in  the  VB  method  the  unknown  parameters  arc  modeled  as  random  variables  and  the  Bayesian  ap¬ 
proach  is  used  to  determine  their  approximate  posterior  densities.  Second,  in  the  EM  method,  each 
M-step  only  receives  the  updated  sufficient  statistics  from  the  E-step,  where  as  in  VB,  as  a  result  of 
the  Bayesian  principle,  the  VBM  step  utilizes  the  updated  priors  of  the  randomized  parameters. 

The  parameters  of  the  clutter  model  and  target  model  are  chosen  as  follows:  Nakagami  parame¬ 
ter  a  =  3,  clutter  power  a2  =  5dB,the  complex  target  response  x  =  1.6  + l.Oi,  the  spectral  density 
of  the  waveform  is  normalized  <!>/  [  =  1,  the  number  of  observations  N  varies  from  10  to  1000. 

Figs.  1.1(a)  shows  that,  for  the  parameters  cr2,  the  variational  algorithm  VBL-PF  outperforms 
the  PXEM  algorithm  when  the  number  of  observations  N  <  200  while  the  VBN-PF  method  outper¬ 
forms  the  PXEM  for  all  values  of  N.  For  the  estimation  of  Nakagami  parameter  a,  the  performance 
of  VB  methods  is  much  better  than  PXEM  in  terms  of  MSE,  as  shown  in  Fig.  1.1(b).  This  is  be¬ 
cause  in  the  PXEM  method,  the  ML  solution  is  obtained  by  solving  a  nonlinear  equation  due  to  the 
non-Gaussian  clutter  model.  For  the  estimation  of  radar  target  response  x,  the  three  methods  have 
very  similar  performance  for  all  values  of  N  as  depicted  in  Fig.  1. 1(c).  The  explanation  for  the  sim¬ 
ilarity  is  that  the  unknown  parameter  is  a  linear  function  of  the  observations  and  follows  a  Gaussian 
distribution.  Note  that  ML  and  Bayesian  algorithms  usually  have  similar  performance  for  Gaussian 
lineal-  models.  The  advantage  of  VB  approach  can  be  observed  for  non-Gaussian  nonlinear  models 
such  as  the  ones  that  involve  a2  and  a.  This  observation  is  consistent  with  existing  literature  on 
variational  Bayesian  studies  [Smidl  and  Quinn,  2008,  Tzikas  et  ah,  2008].  Finally,  The  estimation 
performance  of  the  adjustment  parameter  A  is  not  discussed  as  it  is  only  a  nuisance  parameter  for 
correcting  the  multiplicative  error. 
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Figure  1.1:  Comparison  of  mean  squared  error  (MSE)  between  PXEM,  VBN-PF  and  VBL-PF.  (a) 
clutter  noise  variance  a2,  (b)  Nakagami  parameter  a,  and  (c)  target  response  x. 

1.5  Conclusion 


We  develop  variational  Bayesian  algorithms  for  estimating  multi-parameters  of  a  compound  Gaus¬ 
sian  clutter  model  and  target  response.  The  VB  method  yields  closed  form  expressions  for  the 
posterior  probability  density  functions  and  results  in  improved  estimation  performance  for  the  clut¬ 
ter  model,  especially  for  parameters  of  non-Gaussian  nonlinear  models  and  when  the  number  of 
measurements  is  small. 


Chapter  2 


Range  and  Doppler  Estimation  by 
Weighted  OFDM  Modulation 


2.1  Introduction 

An  increasing  array  of  technologies  in  communications  services  has  resulted  in  scarcity  of  the  fre¬ 
quency  spectrum.  Spectrum  congestion  becomes  a  growing  problem,  thus  limiting  the  operational 
capabilities  of  competing  wireless  systems  due  to  mutual  interference.  The  call  for  spectrum  sharing 
is  to  allow  co-existence  between  various  RF  systems  including  Wi-Fi,  cellular  system,  and  S-band 
radar  that  operate  in  the  same  or  close  frequency  band  while  achieving  normal  operational  perfor¬ 
mance  levels.  Traditionally,  fixed  spectrum  allocation  has  been  used  to  prevent  interference  for 
radar  and  communications  systems  that  operate  in  close  ranges.  However,  with  the  huge  success  in 
wireless  connectivity  and  cellular  communication  system,  more  frequency  bands  are  allocated  for 
commercial  applications,  which  makes  the  fixed  spectrum  allocation  approach  increasingly  difficult 
to  implement. 

One  of  the  alternative  approaches  to  achieve  spectrum  sharing  is  to  design  reconfigurable  wire¬ 
less  platforms,  for  example,  a  single  RF  platform  that  can  be  used  to  realize  communication  and 
radar  functions.  However,  the  intended  radar  functions  are  quite  different  than  communications  sys¬ 
tems.  Radar  is  an  active  sensing  device  that  relies  on  transmission  waveforms  for  target  detection, 
estimation,  and  tracking  with  the  goal  of  achieving  high  resolution.  For  a  communications  system, 
information  bits  transfer  is  accomplished  by  waveform  modulation  and  demodulation.  The  primary 
goal  for  communications  is  spectral  efficiency  and  throughput.  Hence,  one  of  the  challenges  that 
we  face  is  to  design  modulation  waveforms  that  can  be  dual-used  for  radar  and  communication. 
Many  different  modulation  schemes  that  are  traditionally  used  for  communications  can  also  be  uti¬ 
lized  for  radar.  For  example,  spread  spectrum  techniques  have  been  used  in  a  vehicle-to-vehicle 
communication  and  ranging  system  under  a  single  RF  platform  [Uchida  et  al.,  1994].  Orthogonal 
frequency  division  multiplexing  (OFDM)  or  its  variations  [Sit  et  al.,  2011]  have  also  been  used  in 
radar  [Braun  et  al.,  2010], 

There  are  two  noticeable  advantages  of  using  OFDM  modulation.  First,  by  reusing  the  OFDM 
waveform,  limited  extra  hardware  is  needed  for  optimizing  both  wireless  communication  and  radar 
detection.  Second,  recent  research  has  shown  that  OFDM  offers  robustness  performance  against 
multipath  fading  [Sen  and  Nehorai,  2011,  Garmatyuk  et  al.,  2007]  in  highly  mobile  environments 
[Schmidl  and  Cox,  1997].  Note  that  constant  envelope  signaling  with  phase  modulation  is  gener- 
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ally  preferred  for  hardware  realization  [Testi  et  al.,  2013].  However,  the  constant  envelope  OFDM 
waveform  (CE-OFDM)  has  a  relatively  high  peak  to  average  power  ratio  (PAPR).  The  upper  bound 
equals  to  the  number  of  subcarrier  frequencies  [Ochiai  and  Imai,  2001].  Hence,  the  estimation  per¬ 
formance  is  limited  by  the  number  of  independent  individual  OFDM  symbols. 

In  this  work,  we  present  a  non-linear  least  squares  approach  to  estimate  the  velocity  and  range  of 
a  radar  target  based  on  a  weighted  OFDM  (WOFMD)  modulation  scheme.  The  proposed  WOFDM 
method  is  a  non-constant  envelop  modulation  method.  The  weights  of  the  WOFDM  symbols  arc 
designed  by  optimizing  the  error  bound  of  the  estimator  subject  to  the  constraints  of  PARP  and  a 
total  transmission  energy.  Numerical  simulations  show  that  the  weighted  OFDM  scheme  improves 
the  Cramer-Rao  lower  bound  on  the  range  estimation  accuracy  while  controlling  the  maximum  level 
of  PAPR.  Hence,  this  method  provides  a  flexible  and  reconfigurable  realization  mechanism  for  the 
radar/comm  system  co-design  with  an  improved  estimation  range  accuracy  of  radar  targets. 


2.2  Signal  Models 


We  consider  an  OFDM  signaling  system  with  N  subcarriers,  each  modulated  with  a  data  symbol. 
The  time  domain  transmitted  (complex  envelop)  signal  in  the  m-th  pulse  (m  =  0,  •  •  •  ,  M  —  1)  can 
be  written  as 

TV— 1 

sm(t )  =  ^2  aminej27r^ntlm(f,  TPRi,Tp)  (2.1) 

n=0 

where  the  m-th  rectangular  pulse  function  is  defined  as 


1 m{t,  TpRI,  Tp) 


1  m.TpRi  <t<  mTpRi  +  Tp 
0  mTpRj  +  Tp  <  t  <  (m  +  1)Tpri 


(2.2) 


and  am.r,  denotes  the  complex  weights  transmitted  over  the  n-th  subcarrier  for  the  m-th  modulation 
symbol.  For  a  constant  envelope  OFDM  modulation  method,  the  following  condition 


W— 1 

|am,n|2  =  a  =  const.  (2.3) 

71=0 

is  imposed.  Tpri  is  the  pulse  repetition  interval  (PRI).  Tp  is  the  OFDM  symbol  duration  and  the 
subcarrier  spacing  A /  =  .j  =  B/N,  where  B  is  the  signal  bandwidth.  The  individual  subcarrier 
frequency  is 

/„=(/c-f)  +nA/,„  =  °,l,..iV-1  (2.4) 

where  fc  is  the  carrier  frequency  of  the  radar.  In  the  presence  of  a  reflecting  target  at  the  distance 
R  with  the  relative  velocity  of  v,  the  relative  Doppler  shift  is  j3  =  Hence,  the  induced  Doppler 
frequency  at  the  n-th  subcarrier  is 

fn  =  fnP  =  (jc  +  y)2^  (2'5) 

We  assume  that  the  fc  2>  B,  i.e.,  the  carrier  frequency  is  much  larger  than  the  signal  bandwidth, 
hence,  it  is  safe  to  assume  that  the  Doppler  frequency  is  fixed  with  respect  to  the  subcarriers,  i.e., 
the  average  Doppler  frequency  is 

/d  ~  fc—  (2.6) 

c 


10 


thus  applying  the  doppler  shift  on  frequency  fn  — >•  fn  +  f,i  and  the  time  delay  t  — >•  t  —  r  we  have 
the  received  signal 

N—l 

Xm(t )  =  am,neJ'27r{/”+/d){t-T)lm(i,  Iprj,  Tp)  +  wm(t)  (2.7) 

n=0 

We  rewrite  (2.7)  as 

xm(f)  =  H  sm(t  -  T)ej27rfd^T)  +  wm(t)  (2.8) 

where  r  is  the  round  trip  time  delay  between  the  radar  and  target.  A  is  the  target  response.  Note 
that  in  general,  the  reflection  coefficient  A  is  complex  and  may  vary  in  different  subcarrier  channels. 
For  simplicity,  A  is  assumed  to  be  a  deterministic  constant  in  this  work.  wrn  (t)  is  the  clutter  and 
measurement  noise,  following  a  zero  mean  complex  Gaussian  distribution.  Inserting  (2.1),  (2.4)  and 
(2.6)  to  (2.8),  we  obtain 

Xm(t)  =  xsm(t)lm(t,  TPRi,  Tp)  +  wm(t)  (2.9) 

where  the  signal  component  is 

N- 1 

TS  /. \  _  A  SY  j2n(nAf(t-T)-fcT)  j2nfd{t-T)  j2nfct 

Jym\L  )  —  ^  /  j  Ujm,nt:'  c  c 

n= 0 

The  corresponding  baseband  (complex  envelope)  signal  is 

ym(t)  =  xsm(t)lm(t}TpRi,Tp)e~j2nfct  +  wm(t)  (2.10) 

JV— t 

=  Y,  ^«m,ne^(nA/(t-r)-/cT)e^(t-T)lm(f,T Wi,Tv)  +  wm(t)  (2.11) 

n= 0 

Note  that,  for  an  OFDM  system,  the  sampling  interval  Ts  is  typically  chosen  such  that  A fTs  = 
1/N.  Hence,  the  sampled  version  of  ym(l )  with  (Z  =  0,  •  •  •  ,  N  —  1)  becomes 

N—l 

ym (Z)  =  A  Y  «m/2^!e“i2"(nA/+/c)T  x  ej2nfd ^s{l+mQN)-r)  +  Wm^  (2.12) 

n=0 


where  Q  =  TRRi/Tp.  Next,  we  define 

Ys  =  Aam  ne_727r(nA/+/c)rei27r/d(TsmQiV)  ^  13) 

Using  the  definition  of  inverse  DFT  and  Ts  =  XYJ'  we  can  re_wr'te  (2-12)  as 

ym{l)  =  (V{IDFT  [Y^n]  }ej2n^T^e-^dT  +  wm(l )  (2.14) 

o'27T  f  d  — 

the  term  amounts  to  modulation  in  the  time  domain.  Hence  the  DFT  of  ym(l)  has  a 

corresponding  frequency  shift  Thus,  applying  DFT  on  (2.14)  implies 


Y'  =  NYS 

m.n  * 


e~j2nfdr  +  ^ 


(2.15) 


In  radar  estimation,  we  can  assume  that  f,i  <  A /  thus  the  only  term  that  is  influenced  by  modulation 

is  e-727r(nA/+/c)r  wjq1  n  —  n  _  JjL  and 


^rn,n  —  A" &m,n& 


-jM(n-U)Af+fc)Tj2nfd(TsmQN)  v  p-j2nfdr 


(2.16) 
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which  can  be  simplified  as 


Y^n  =  +  Wm,n  (2.17) 

where  A'  =  ANe~i‘2'K*cT  and  is  invariant  term  respect  to  m  and  n.  Wrnjl  is  the  DFT  of  the  noise 
term  wm(l).  Dividing  the  known  OFDM  symbols  arn,n  on  both  sides  of  (2.17),  we  obtain, 

=  A'ejMmfdTsQN-nAfr)  +  ^  (2.18) 


where  Y,„  „  = 


YL 


and  W° 


wm,„ 

O' m,n 


2.3  Minimum  Nonlinear  Least-Squares  Estimator 

We  let  A'  =  bej,IJA  be  the  target  response  with  b  and  (1>a  being  the  magnitude  and  phase  of  target 
response.  The  parameter  vector  of  interest  is 

0  =  [b,  c/)A,  fd,  t\  (2.19) 


From  the  signal  model  (2.18),  the  nonlinear  least-squares  error  function  is, 

M—l  N—l 


Ln(Y-,e )  =  ]T  J2\W™,n\2 
m= 0  n=0 
M—l  N—l 

=  ££  | Ym  n  —  be3<l>Aej2n(mfdTaQN-nAfT)\2 


m= 0  n=0 

where  M  is  the  number  of  pulses  in  the  coherent  processing  interval.  Hence,  we  obtain 

M—l  N-l 

Ln(Y ;  6)  =  J2  (\Ym,n\2  +  b2  -2M{Y*hnbej*Aej2w(mfdTsQN-nAfT) 

771=0  71=0 

The  par  ameters  6  can  be  determined  by  minimizing  the  error  function,  i.e., 

find  6  =  argmin  Ln(Y;9 ) 

0 


(2.20) 


(2.21) 


(2.22) 


Since  there  are  four  parameters  in  the  vector  6 ,  a  sequence  of  optimization  steps  arc  taken.  We  staid 
by  taking  partial  derivative  of  the  error  function  with  respect  to  the  unknown  phase  4>a-  which  yields 


dLn(Y-6) 

d<P  A 


M- 1  N-l 

=  2b^2Y^  (23fH  -  jYm,nbe-^Ae-j2^mfdTsQN-nAfT)])  =  0  (2.23) 

771=0  71=0 


Note  that  the  2-D  discrete  time  Fourier  transform  of  Ym  ,,  is 


M—l  N-l 

Z(fd,r )  4,I,eitaA/T)e‘j2m/jTsQiV 

771=0  71=0 

Inserting  (2.24)  in  (2.23)  we  obtain 

X{-jbe-MAZ(fd,T)}  =  0 


(2.24) 


(2.25) 
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Since  b  and  \Z(fd,  t)|  are  positive  real  values,  only  the  real  component  of  the  phase  term  must  be 
zero.  Hence  we  have 

cos  (i  +  ~  r')) =  0  (-2‘26') 

This  is  a  simple  trigonometric  equation  and  the  solution  is 

4>a  =  kir  +  ZZ(fd,r),  (2.27) 


where  k  is  an  integer.  Using  this  estimate  (2.27)  in  the  error  function  (2.21)  and  ignoring  the  terms 
independent  of  b,  we  obtain 

Ln(Y ;  (f>A,  b,  fd,r)  =  2b‘R{e~zz^T)  Z(fd,r)}  -  MNb2  (2.28) 


which  simplifies  to 

Ln(Y ;  fa,  b ,  fd,  t)  =  2b\Z(fd,  r)|  -  MNb 2 

By  equating  the  partial  derivative  with  respect  to  b  to  zero  yields 

i  \Z(fd,T)\ 

MN 


(2.29) 


(2.30) 


Inserting  the  estimate  (2.30)  in  (2.29)  we  obtain 

Ln(Y;  <j>A,  b,  fd,r)  =  (2-31) 

which  is  a  scaled  the  periodogram  of  the  2D-DTFT  in  (2.24).  Hence,  we  obtain  the  final  optimiza¬ 
tion  problem  for  unknown  delay  and  doppler  frequency  parameters  as 

find  fd,T  =  argmax  (2.32) 

fd,r  MN  v  ' 

The  solution  to  (2.32)  can  be  determined  by  finding  the  peak  of  un-windowed  2D  periodogram  of  the 
signal  Fm>n.  This  estimation  can  be  accomplished  by  a  two-dimensional  search  [Rife  and  Boorstyn,  1974]. 


2.4  Cramer- Rao  Lower  Bounds 


Note  that  if  Wrn,n  in  (2.18)  is  a  Gaussian  random  process,  0  is  a  large  sample  realization  of 
the  maximum-likelihood  estimator  of  6.  Since  the  maximum-likelihood  estimator  is  expected  to 
achieve  the  Cramer-Rao  Low  Bound  (CRLB)  as  M  or  N  increases,  it  follows  that  under  the  Gaus¬ 
sian  assumption  no  other  estimator  could  perform  better  in  large  samples  than  6.  In  this  section,  we 
derive  the  Cramer-Rao  lower  bound  of  6.  The  probability  model  of  the  measurements  Y  =  { Ym/n  } 
is 


M— 1  N—l 


=  11 II 


m= 0 n— 0 


7 raf. 


-exp 


*"m, n\ 


I  Ymn-  bej*Aej 2^(™/dTsQiV-nA/r)|2  |2  j  (2_33) 


The  log-likelihood  function  is 


M—l  N—l 

l(Y\9)  =  +  b2\am,n\2  -  23F i{Y^nbej^Aam,ne^mfdT^N-nA^})  (2.34) 

m=0 n= 0 
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The  Cramer-Rao  bounds  for  the  parameter  estimates  arc  given  by  the  diagonal  elements  of  inverse 

of  the  4x4  Fisher  information  matrix  J  [Kay,  1993], 

Var(0i)  >  [J-1];;  (2.35) 


where  8,  arc  the  four  scalar  parameters  in  the  vector  6  in  (2.19)  with  i  =  1,  2, 3, 4. 

Due  to  space  limitations,  a  detailed  derivation  of  the  elements  arc  omitted  in  this  work.  We  only 
present  the  final  results  here.  The  determinant  of  Fisher  information  matrix  J  is 

256  56 

det(J)  =  — =-( MaTT2TsQNAf)2(MQM  ~  P&)(aQN  -  P'n)  (2-36) 

(T° 

where  PM  =  £^=o  m  =  M (M  -  l)/2,  QM  =  £^=o  =  M(M  ~  l)(2M  ~  l)/6,  Qn  = 

J2n=o  r>‘2\am-n\2  and  PN  =  ri\am/n\2.  The  Cramer-Rao  low  bound  on  the  estimate  of  b  is 

the  first  element  along  the  diagonal  of  J-1  given  by 

1  128b6 

crb6  =  det(J)  (MQm  ~  Pm)(Qn  ~  pn)—^- 

Man4{TsQN)2{Af)2  (2.37) 


The  diagonal  elements  of  the  J  1  are  evaluated  using  the  adjoint  method.  Next,  using  (2.36)  and 

2 

upon  simplification  we  obtain  CRH/,  =  jp^h-  The  CRLB  on  4>a  is  the  second  element  along  the 
diagonal  of  J-1  and  is  obtained  as, 


CRB  =  ^  MciQmQn  -  Pmpn 

^  b 2  2Mcl(MQm  ~  P2M){aQN  -  P2) 

The  CRLB  on  f,i  is  the  third  element  along  the  diagonal  of  J-1  given  by, 

CRB^  =  dit(J)f: ^M)\aQN-pi) 

Again  upon  simplification  we  obtain 

rRR  =  _ _ 

fd  8 abH2(TsQNY{MQM  ~  Pi) 


(2.38) 


(2.39) 


(2.40) 


Finally,  the  CRLB  on  r  is  the  forth  element  along  the  diagonal  of  J  1  defined  as, 

CRB^  =  a^%KJ^aPt‘qn?(MQm -  Pll) 

Hence,  upon  simplification  we  obtain 


CRBr 


ac r2 

2 (aQN  -  P2  )62M(27rA/)2 


(2.41) 


(2.42) 
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2.5  Weighted  OFDM  Symbol  Design 


We  use  the  classic  Cramer-Rao  lower  bound  (CRLB)  to  evaluate  the  performance  of  the  nonlinear 

least-squares  estimation.  Eqn.  (2.42)  reveals  that  the  CRLB  on  time  delay  estimate  t  is  determined 

2 

by  the  factor  W,  a  term  related  to  SNR,  the  number  of  pulses  M,  the  frequency  separation  A / 
and  the  factor  —  -  a  2  which  depends  on  OFDM  symbol  amplitudes  I am  n\.  Similarly,  Eqn.  (2.40) 

shows  that  CRLB  on  the  Doppler  frequency  estimate  f(i  depends  on  the  SNR  value  ,  the  constant 
M  2  and  the  sampled  length  of  the  pulse  repetition  interval  TSQN,  but  independent  of  the 

1V1 LJ  m  -*  m 

spectral  shape  of  OFDM  symbols  |am,n|- 

By  the  definition  of  CRLB,  we  know  that  the  determinant  of  the  Fisher  information  matrix 
det(J)  in  (2.36)  is  the  common  denominator  for  each  of  the  CRLBs  given  in  (2.37),  (2.38),  (2.40) 
and  (2.42).  The  aim  for  weighted  OFDM  is  to  design  OFDM  symbols  to  minimize  the  inverse  of 
det(J)  subject  to  proper  constraints,  i.e., 

find  am,n  =  arg  min  1  ,  (2.43) 

am,n  ClLj —  Tyy 

{(Cl)  min(|am,in|)  >  r/0  max(|amin|) 

(C2)  PAPR  <  r/i  (2.44) 

(C3)  £n=o  |am,n|2  <  a 


Condition  (Cl)  is  to  impose  a  constraint  on  the  lower  bound  of  the  cost  function  aQ^_P2  by 


min(|amn|) 

- 77 - TT  <  VO  <  1 

max(|amj„|) 


(2.45) 


The  purpose  of  this  constraint  is  to  avoid  bandwidth  loss.  Otherwise,  OFDM  symbols  am,n  would 
allocate  zero  power  in  most  of  sub-bands  n,  causing  bandwidth  loss  in  the  transmission  signal.  As 
a  result,  the  range  resolution  could  be  severely  degraded.  Condition  (C2)  is  to  set  an  upper  bound 
r/i  <C  N  on  the  peak  to  average  power  ratio  (PAPR)  defined  by 


PAPR  =  a  1  max  |sm(t)|2 

0  <t<Tp 


(2.46) 


because  a  high  PAPR  imposes  severe  burden  on  the  transmitter  due  to  limited  amplification  range  of 
the  RF  amplifier.  For  instance,  when  phase  shift  keying  (PSK)  is  used,  the  upper  bound  on  PAPR  is 
N  [Ochiai  and  Imai,  2001],  Finally,  Condition  (C3)  is  the  total  power  constraint.  The  optimization 
problem  defined  in  (2.43)  and  (2.44)  is  solved  numerically  using  convex  optimization  by  the  active- 
set  constrained  nonlinear  method  [Han,  1977]. 


2.6  Simulation  Results 

The  estimation  performance  of  the  proposed  weighted  OFDM  (WOFDM)  is  compared  with  the 
QPSK  based  constant  envelope  OFDM  (CE-OFDM)  scheme.  For  QPSK,  we  choose  |amjn|2  = 
1/N.  The  phase  of  the  OFDM  symbols  am^n  is  a  uniformly  distributed  random  variable.  Under 
this  modulation  scheme,  the  cost  function  clQn  —  P ^  =  jV^2~1  is  a  constant.  The  CRLB  of  the 
Doppler  in  (2.40)  is  independent  of  modulation.  The  CRLB  of  delay  estimate  for  CE-OFDM  is  also 
independent  of  symbol  amplitude  |amjn|  and  is  given  by 

CRBt  =  3a2/  (2 (N2  -  1)62M(ttA/)2)  (2.47) 
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For  WOFDM,  both  amplitude  and  phase  of  the  weights  aTO)„  arc  derived  from  the  optimization 
problem  (2.43)  and  (2.44).  In  the  simulations,  the  bandwidth  B  =  10  MHz,  the  number  of  sub¬ 
carriers  is  N  =  128,  and  the  carrier  frequency  is  fc  =  2  GHz.  The  frequency  separation  is  A /  = 
781.25  kHz.  The  signal  energy  is  normalized  (i.e.,  a  =  1)  and  the  SNR  varies  from  1  ~  13  dB.  The 
pulse  repetition  interval  is  Tppj  =  NTp.  The  resolution  of  the  Doppler  spectrum  is  1  / (TSQN M). 
The  delay  spectrum  has  a  resolution  of  Tppi/N  =  Tp. 

The  true  values  of  parameters  of  a  simulated  moving  target  arc  given  as  follows.  The  Doppler 
frequency  is  fcj  =  ,5/Tppj ,  time  delay  is  r  =  Tppj/2  and  the  number  of  pulses  is  M  =  50.  The 
optimization  parameters  are  set  as  //q  =  0.2  and  r/i  =  3.  The  estimates  }'cj  and  t  arc  obtained  by  the 
nonlinear  least-squares  estimator  formulation  in  (2.32).  By  the  active-set  constraint  optimization 
method,  we  calculate  the  weights  am>n  numerically  [Han,  1977].  Table  2.1  presents  the  statistics  of 


Table  2.1: 

Variations  of  PAPR  Values 

Modulation 

CE-OFDM 

WOFDM 

average  PAPR 

5.643 

3.004 

variance 

1.348 

0.016 

PAPR  from  the  CE-OFDM  vs.  WOFDM.  The  WOFDM  method  yields  a  more  stable  (i.e.,  smaller 
variance)  and  lower  PAPR  value  compared  to  the  CE-OFDM  method. 

Fig.  3. 1  depicts  the  theoretical  CRLB  and  the  numerical  values  for  the  delay  and  Doppler  esti¬ 
mates.  For  the  delay  estimate  in  Fig.  3.1(a),  the  plots  show  that  WOFDM  modulation  given  in  (2.42) 
achieves  a  smaller  CRLB  than  CE-OFDM  given  in  (2.47).  This  is  because  by  varying  the  magni¬ 
tudes  of  the  OFDM  symbols,  we  reduce  the  value  of  the  max(|sm(f)[),  i.e.,  bringing  it  closer  to  the 
mean  value  of  (|sm(f)|).  Hence  we  achieve  a  higher  value  for  det(J),  thus  leading  to  a  lower  CRLB 
for  delay  estimate.  As  expected,  for  Doppler  estimate,  the  CRLBs  arc  the  same  for  the  two  mod¬ 
ulation  schemes.  Note  that  the  variance  plots  of  the  parameter  estimates  by  the  WOFDM  method 
closely  follows  the  variance  plots  of  CE-OFDM  only  after  the  SNR  goes  above  8  dB.  This  phe¬ 
nomenon  is  often  called  threshold  effect  in  nonlinear  estimation  problem  [Rife  and  Boorstyn,  1974]. 
It  means  at  low  SNR,  there  is  usually  a  range  of  SNR  in  which  the  mean-squared  error  (MSE)  rises 
rapidly  as  SNR  decreases.  In  other  words,  only  at  high  SNR,  the  computed  SNR  follow  the  CRLB 
closely,  as  shown  in  Fig.  3.1.  Fig.  3.2  represents  the  improvements  in  CRLB  using  the  WOFDM 
compared  to  the  CE-OFDM  with  respect  to  different  values  of  r/o  and  r)\.  The  plots  show  the  trade¬ 
off  between  the  CRB,  PAPR  and  bandwidth  preservation,  i.e.,  it  is  not  possible  to  achieve  low  CRB 
and  low  PAPR  simultaneously. 

2.7  Conclusion 

A  weighted  OFDM  scheme  is  proposed  which  provides  a  lower  PAPR  compared  to  the  QPSK  based 
constant  envelope  OFDM  scheme.  Theoretical  analysis  and  numerical  simulation  demonstrate  that 
the  CRLB  bound  for  Doppler  (velocity)  estimate  does  not  change  with  modulation,  however,  the 
CRLB  of  delay  (range)  estimate  is  improved  with  the  weighted  symbols  by  WOFDM.  The  pro¬ 
posed  WOFDM  scheme  provides  a  promising  means  to  achieve  co-existence  between  radar  and 
communications  via  a  reconfigurable  RF  platform. 
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Delay  (tau):  variance  comparison  Doppler  (fd):  variance  comparison 


(a) 


(b) 


Figure  2.1:  Computed  variance  and  CRLB  for  CE-OFDM  vs.  WOFDM  modulation,  (a)  Delay 
estimate  f .  (b)  Doppler  frequency  estimate  fci 


Figure  2.2:  CRLB  for  WOFDM,  a  compromise  between  the  upper  bound  (CE-OFDM),  lower  bound 
and  PAPR 
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Chapter  3 


Waveform  Encoding  for  Nonlinear 
Electromagnetic  Tomographic  Imaging 

3.1  Introduction 

Electromagnetic  (EM)  tomographic  imaging  is  an  inverse  scattering  problem  which  has  a  wide 
range  of  applications  in  medical  imaging,  geophysical  exploration,  nondestructive  testing,  and  tar¬ 
get  identifications  [Chew  and  Wang,  1990,  Franchois  and  Pichot,  1997,  Isernia  et  ah,  1997].  In  EM 
tomographic  imaging,  source  antennas  transmit  EM  signals  into  a  medium  under  test  and  receive 
scattering  signals.  Based  on  the  underlying  Maxwell’s  equations,  inversion  methods  arc  employed 
to  reconstruct  a  spatial  distribution  of  material  parameters  such  as  the  dielectric  permittivity  and/or 
magnetic  permeability  of  the  target  and  the  surrounding  medium,  thus  turning  recorded  scattering 
data  into  images  [Fhager  et  al.,  2006,  Liu  et  al.,  2002]. 

Mathematically,  electromagnetic  tomographic  imaging  is  an  inverse  problem,  which  can  be 
written  as 

Vj  =  -Aj(p( r);  Sj)  +  r/J;  j  =  1,  ■  ■  ■  ,  Q  (3.1) 

The  goal  is  to  infer  model  parameters  (i.e.,  material  property  values  p(r),r  €  S2)  from  the  ob¬ 
served  data  (i.e.,  measurement  y]  €  Oil  x  [0,  T])  in  response  to  the  j-th  excitation  source  Sj  based 
upon  a  Maxwell  wave  model,  fl  denotes  the  imaging  field,  T  is  the  integration  time  of  the  signals 
measured  at  the  receivers  positioned  at  the  boundary  Oil  of  the  imaging  field.  Q  is  the  number 
of  excitation  sources.  A:) ( • )  is  the  nonlinear  operator  determined  by  the  wave  model  and  rjj  is  the 
noise  or  disturbance  term.  The  image  to  be  reconstructed  is  a  spatial  distribution  of  the  parame¬ 
ter  p.  The  dimension  of  p,  determined  by  the  discretization  of  the  image  field  fl,  is  often  much 
larger  than  the  dimension  di > |  (i.e.,  the  number  of  receivers)  for  data  y} .  Hence,  the  full-wave 
inverse  problem  (3.1)  is  considered  ill-posed,  often  requiring  regularization  of  the  inverse  problem. 
For  example,  the  well  known  compressive  sensing  method  would  impose  a  sparsity  constraint  on 
p  [Candes  and  Tao,  2006,  Lustig  et  al.,  2008],  thus  effectively  reducing  the  dimensionality  of  the 
solution  space  for  p. 

The  imaging  problem  (3.1)  we  consider  is  more  challenging  in  that  it  involves  solving  a  non¬ 
linear  full  wave  inverse  problem  with  active  excitation  sources.  Classic  approaches  to  solving  (3.1) 
include  the  least  squares  optimization  method  [He  et  al.,  1997,  Jia  et  al.,  2002]  or  the  iterative  New¬ 
ton’s  method  [Dorn  et  ah,  1999].  In  both  cases,  iterative  computational  techniques  are  employed. 

The  cost  of  computation  depends  on  the  size  of  the  data  volumes  (i.e.,  the  number  of  transmis¬ 
sion  sources  and  the  number  of  receivers)  and  on  the  discretization  of  the  wave  model.  The  com¬ 
plexity  of  computing  gradient  and  Newton  updates  -  aside  from  issues  with  non-convergence  and 
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non-uniqueness  -  become  the  major  impediment  withstanding  successful  application  of  full  wave  in¬ 
version  to  industry-size  data  volumes,  for  example,  large  scale  seismic  imaging  [  Krebs  et  al.,  2009] 
and  high  resolution  medical  imaging  [Fhager  et  al.,  2006]. 

In  this  work,  we  will  address  this  issue  by  means  of  dimensionality  reduction  through  multi¬ 
ple  source  waveform  encoding.  Typical  full  wave  tomographic  imaging  operates  in  a  single-input 
multiple-output  (SIMO)  configuration.  An  image  is  reconstructed  from  measured  data  in  response 
to  a  single  excitation  antenna  source.  The  reconstruction  process  continues  till  all  the  sources  arc 
excited.  For  large-scale  imaging  such  as  seismic  exploration,  the  number  of  EM  sources  is  very 
large.  Not  only  the  computational  cost  is  high,  the  operational  expenditures  of  each  data  collection 
process  is  also  significant.  In  this  case,  multiple  source  excitation  becomes  appealing.  For  exam¬ 
ple,  in  seismic  imaging,  multiple  sources  arc  excited  simultaneously  to  form  a  so-called  supershot 
to  probe  the  imaging  field  [Godwin  and  Sava,  ],  which  means  the  imaging  configuration  becomes 
multiple-input  multiple-output  (MIMO). 

The  recorded  data  arc  processed  to  form  an  image.  The  image  is  updated  when  new  measure¬ 
ment  data  is  available.  This  procedure  is  repeated  until  the  image  converges  or  a  pre-determined 
stopping  criterion  is  met.  However,  multiple  wave  simultaneous  excitation  induces  cross-talk  noise 
due  to  wave  interference,  which,  if  not  treated,  will  cause  image  artifacts.  Therefore,  signal  pro¬ 
cessing  techniques  such  as  waveform  encoding  arc  needed  to  mitigate  cross-talk  noise  in  the  image 
in  order  to  achieve  high  quality  imaging  while  reducing  the  computational  complexity. 

The  contribution  of  this  work  is  threefold.  First,  we  develop  three  different  waveform  encoding 
techniques,  i.e.,  random  phase  encoding,  waveform  delay  encoding,  and  uniform  weight  encoding 
for  the  full-wave  EM  imaging  problem.  We  show  that  the  random  phase  encoding  method  results 
in  constant-envelop  waveforms  and  produces  the  best  performance  in  terms  of  convergence.  Second, 
this  work  extends  our  early  work  on  microwave  imaging  in  a  SIMO  configuration  [Dong  et  al.,  2012]. 
We  show  that  using  simultaneous  sources  made  of  superposition  of  encoded  sources  is  able  to  ac¬ 
celerate  iterative  algorithms  for  electromagnetic  full-wave  inversion,  thus  demonstrating  the  effec¬ 
tiveness  of  waveform  encoding,  a  common  signal  processing  technique,  to  improve  computational 
efficiency  of  classic  nonlinear  inverse  problems.  Third,  although  waveform  encoding  techniques 
have  been  studied  for  acoustical  wave  imaging  [  Krebs  et  al.,  2009]  and  our  early  work  on  MIMO 
ultrasonic  imaging  [Dong  and  Jin,  2013,  Dong  et  al.,  2014],  there  is  still  a  lack  of  research  for  EM 
full-wave  imaging  in  applications  where  the  use  of  EM  wave  energy  is  critical.  In  this  work,  we  de¬ 
velop  iterative  algorithms  that  solve  time  domain  Maxwell’s  equations  with  coupled  electric  fields 
and  magnetic  fields  using  waveform  encoded  excitations  and  demonstrate  faster  convergence  com¬ 
pared  with  our  prior  SIMO  imaging  algorithm  in  [Dong  et  al.,  2012]. 

3.2  MIMO  Imaging  Problem  Formulation 

Consider  the  time  dependent  Maxwell’s  equations  in  an  isotropic  medium  [Chew,  1994] 

e(r)  dEQt'  ^  -  v  x  H(r,t)  +  cr(r)E(r,t)  =  -J(r,t )  (3.2) 

Mr)  +  v  x  E(r,t)  =  0  (3.3) 

where  the  quantities  E( r,  t)  =  (Ex,  Ey.  Ez,  t )  and  H( r,  t)  =  ( Hx ,  Hy .  Hz.  t)  are  the  electric  and 
magnetic  field  intensities  defined  on  O  x  [0,  T],  where  r  S  f!is  the  imaging  region.  The  material 
parameters  including  electrical  permittivity  e(r),  magnetic  permeability  /./(r)  and  conductivity  cr(r) 
are  independent  of  time.  The  quantity  J(r,  t)  =  ( Jx .  Jy .  Jz,  t )  is  the  electric  current  density  of  an 
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Figure  3.1:  2D  MIMO  tomographic  imaging  geometry.  At  the  m-th  transmission,  the  m-th  group 
of  multiple  current  sources  are  excited  to  emit  TM  waves  to  interrogate  the  imaging  field.  Scattered 
signals  are  received  by  all  the  sensors. 


external  charge  and  is  considered  as  the  external  source.  The  initial  conditions  are  E(r,  0)  =  0  and 
H(  r,0)  =  0. 


3.2.1  2D  EM  tomographic  imaging  problem 

We  consider  a  2-D  inverse  scattering  problem  where  only  the  transverse  magnetic  (TM)  waves  are 
studied.  In  this  case,  Ez  ^  0,  Hz  =  0.  Furthermore,  for  the  sake  of  simplicity,  we  assume  that  only 
the  Ez  component  is  recorded,  i.e.,  we  let  Ex  =  Ey  =  0.  By  omitting  the  space  parameter  r  for  the 
moment,  the  Maxwell’s  equations  in  (3.2)  and  (3.3)  become 


dEz 

dt 

dHy 

dt 


( t ) 

(t) 


——— Ez(t )  + 


1  dEz 


H  dx 


(<)> 


and 


dHv  dHr\  .  .  1 

-w+nf)  {t>-i 

dHx  .  ,  1  dEz 


Jz(t ) 

C t ) 


(3.4) 

(3.5) 


where  Jz(t)  is  the  current  source  polarized  in  the  z-dircction  emitted  from  the  transmit  antennas. 
Specifically,  in  the  tomographic  imaging  setup  shown  in  Fig.  3.1,  the  excitation  signal  source  can 
be  expressed  by 

Jz,j(r,t)  =  Sj(t)8(r  -  r*),  J  =  1,  2,  •  •  ■  ,Q  (3.6) 

where  r1-  is  the  position  of  the  j-th  transmit  antenna.  Q  is  the  total  number  of  transmitters.  Sj  (t) 
is  the  excitation  signal  emitted  by  the  j-th  transmitter.  d(r)  is  the  Dirac  delta  function.  The  short 
pulsed  signal  generated  by  a  source  illuminates  the  imaging  field  and  the  unknown  target.  The 
transient  field  data  are  collected  at  the  receiver  positions  Furthermore,  we  assume  in  this  work 
that  the  transmitter  and  the  receiver  positions  are  the  same,  i.e.,  r’j  =  rf-. 

Next  we  combine  Lm  individual  excitation  sources  at  the  m-th  transmission,  leading  to  the  m-th 
“supershot” 

Lm 

Jzj(r*,£),  m  =  1,  •  •  ■  M  (3.7) 

i= 1 

which  is  the  combined  Lm  >  1  sources.  Note  that  the  total  number  of  sources  remains  unchanged, 
i.e.,  Q  =  Ylm=i  Lm-  When  Lm  =  L  is  a  constant,  then  M  =  Q/L.  It  is  easy  to  see  when  Lm  =  1, 
the  MIMO  configuration  reduces  to  the  SIMO  configuration. 


20 


Next,  let  the  symbol  p  =  [e  a  //  7  denote  the  parameter  vector  which  stands  for  the  material 
properties  of  the  internal  structure  of  the  medium  (or  objects)  to  be  reconstructed.  The  imaging 
problem  is  to  seek  the  solution  of  the  nonlinear  equations 

Ez,m(rrj,t)  =  EZtm(p;  rrj,  t)  +  (3.8) 

where  m  =  1,  •  •  •  ,  M  is  the  index  for  transmission  sources  and  j  =  1,  •  •  •  ,  Q  is  the  index  for 
receivers.  L2,mfr’\  t)  is  the  measured  z-componcnt  of  the  electric  field  at  the  receiver  antenna 
position  r'j.  Ez>m(p-,rrj,t)  is  the  calculated  electric  field  for  the  parameter  p  at  the  j-th  receiver 
position  r'j  in  response  to  the  m-th  excitation  source  ./2-r„  (rk  t).  The  calculation  of  Ez/m  (p:  r!j,  t)  is 
based  upon  the  underlying  Maxwell’s  equations  (3.4)  and  (3.5).  Next,  organizing  all  the  measured 
electrical  field  data  into  a  space-time  matrix  with  respect  to  (rr,  t) 

em  =  [Ez,m(ri,t),  Ez,m(r2,  f),  • ' '  ,EZim(rrQ,t)]T  (3.9) 

Accordingly,  the  imaging  problem  (3.8)  can  be  written  as  a  system  of  nonlinear  operator  equations 

Gm  —  AmijP'i  J z.rn  1 1'  ;())  T  rjm  (3.10) 

where  the  j-th  component  [Am(p;  Jz,m.( r*,  t))]j  =  Ezzm  (p:  r'-.  t)  is  the  calculated  electrical  field  at 
the  j-th  receiver  antenna.  rjm  is  the  additive  noise  or  disturbance.  The  problem  of  interest  is  to  find 
out  the  solution  p( r)  for  (3.10)  over  the  imaging  region  r  G  0  given  proper  initial  condition  and 
boundary  conditions  and  waveform  encoded  excitation  signals  J2;m(r*,  t). 

3.2.2  Inversion  algorithm  by  adjoint  fields 

We  use  iterative  Newton’s  method  to  calculate  the  material  property  value  p.  Let  piU>  be  the  initial 
approximation,  the  update  equation  is  given  by 

pk+l  =  pk  +  p8pk  (3.11) 

where  7  is  a  relaxation  factor.  5pk  is  the  increment  value  at  k- th  iteration.  We  employ  the  adjoint 
method  to  calculate  5pk  [Dong  et  ah,  2012],  which  takes  the  form  of 

5pk  =  (^Xm{pk,JZ:m(rt,t))^  d k,m  =  k  mod  M  (3.12) 

where  d(fc)  =  em  -  e$\ anx[0)T]  is  the  difference  signal.  e$\dnx[o ,t]  =  Am{p(k\  Jz,m{*\t))  is 

the  calculated  electric  field  signal  value  at  the  boundary  at  the  k- th  iteration  in  response  to  the  m-th 
transmission.  Here  the  symbol  Am  ( ■ )  stands  for  a  locally  uniformly  bounded  Frechet  derivative  of 
A and  Am{-)*  is  the  adjoint  operator  of  Am{  j  [Dorn  et  ah,  1999,  Dong  and  Jin,  2013].  The 
symbol  mod  stands  for  the  modulo  operation. 

The  outline  of  the  algorithm  is  given  in  Table  3.1.  The  inversion  algorithm  involves  four 
steps.  Step  1  solves  a  forward  problem,  i.e.,  from  A:-th  value  pk,  we  calculate  the  projected  value 
Am{pk]  J2,,n(r*,t))  at  the  boundary  of  the  imaging  field  <9Q,  which  yields  ekn.  Step  2  calcu¬ 
lates  the  difference  signal  dk.  Step  3  is  the  inversion  algorithm  that  calculates  the  adjoint  fields 
( Am(pk ;  Jz,m{ r*,  i)))  -  Step  4  computes  the  update  and  evaluates  the  stopping  criterion.  From  the 
operator  point  of  view,  we  can  show  from  (3.7)  that  Am(p ;  Jz^mi^.t))  =  i.  Aj(p;  J2j(r*-,  t)). 
We  omit  the  proof  due  to  space  limitations. 
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Table  3.1:  Full  Wave  Inversion  Algorithm  by  Adjoint  Fields 

Initialize 

p  <—  p° 

Iteration 

while  \\pk+l  —  pk  2  >  e  do 

1.  ekm  =  Arn(pk) 

//  forward  problem 

2.  dk  =  em  -  ekn 

//  difference  signal 

3.  5pk  =  (A'm(pk))*  dk 

//  adjoint  fields 

4.  pk+1  =  pk  +  yd pk 

//  update 

end 

3.3  Waveform  Encoding  in  EM  Imaging 

MIMO  waveform  signaling  is  well  studied  in  signal  processing  community  to  improve  spatial  diver¬ 
sity  for  better  estimation  and  detection  [Jin  et  al.,  2010].  In  this  work,  the  goal  of  MIMO  signaling 
is  to  improve  computational  efficiency  of  nonlinear  imaging  by  accelerating  the  convergence  of 
the  iterative  algorithm  while  achieving  high  resolution  images.  We  consider  a  Gaussian  modulated 
pulse  as  the  basic  excitation  signal,  which  takes  the  form  of  Sj(t)  =  s(t )  =  cos(27r/c£)e-27r""t~/Tp 
where  tp  is  the  pulse  width.  Next,  we  examine  three  waveform  encoding  techniques  for  imaging. 

3.3.1  Waveform  encoding  design 

(1)  Random  phase  encoding.  The  method  applies  a  random  phase  term  Wj  =  el0j  to  each  excitation 
signal  Sj(t),  where  i  denotes  the  imaginary  unit  and  6j  =  Uniform(0,  27t).  In  particular,  we  consider 
binary  phase  shift  keying  scheme  where  only  { ein  =  —1}  and  {e'°  =  +1}  are  chosen  randomly 
following  a  Bernoulli  distribution.  Hence,  the  aggregated  excitation  signal  takes  the  form  of 

Lm 

t)  =  ^2  WjSjifiSi r*  -  r$),  m  =  1,  ■  •  •  M  (3.13) 

3= 1 

Our  simulation  studies  reveal  that  this  method  reduces  cross-talks  while  achieving  accelerated  con¬ 
vergence. 

(2)  Time-delays  encoding.  To  reduce  the  cross-talks,  a  fixed  time  delay  ry  is  applied  to  the  Lm 
sources  at  the  m-tli  transmission. 

Lm 

Jz%(r\t)  =  ^sj(t-rj)s(rt  ~rtj),  rn  =  l,---M  (3.14) 

3= 1 

where  r:/  is  the  relative  delay  for  the  j-th  source.  A  good  choice  of  time  delay  is  the  pulse  width, 
i.e.,  Tj  =  tp,  for  Gaussian  modulated  pulses. 

(3)  Uniform  weight  encoding.  In  this  case,  Lm  >  1  excitation  signals  are  emitted  simultaneously 
from  Lm  antennas,  the  aggregated  electric  current  density  becomes 

Lm 

=  ~ri)>  m  =  (3.15) 

3= i 


22 


However,  the  mutual  interference  between  the  excitation  signals  may  cause  cross-talks,  which  de¬ 
grades  the  image  quality. 


3.3.2  Impact  of  excitation  sources 


To  examine  the  impact  of  the  excitation  sources  on  the  solution  of  the  imaging  problem,  we  re-write 
the  Maxwell’s  equation  as  a  second  order  wave  equation.  Stalling  from  the  Farday’s  Law  in  the 
vector  form,  Vx£  =  -/j®,  and  take  the  curl  of  both  sides,  we  obtain  VxVxf  =  — ft  J^(V  xH). 
Note  that  V  x  V  x  E  =  —  V2E,  and  using  (3.2),  we  obtain 

^  («*  +  a^)  (3-16) 

Next,  consider  the  TM  mode  where  only  Ez  is  to  be  reconstructed,  we  can  re-write  (3.16)  as  a 
second-order  scalar  wave  equation 


d2Ez  282EZ 
dt 2  C  dz 2 


(3.17) 


where  the  forcing  term  in  (3.17),  f(z,  t)  =  Ijjj  ( Jz  +  aEz),  depends  on  the  excitation  source  Jz. 
More  importantly,  by  the  theory  of  wave  equation  [Evans,  2010],  the  solution  to  inhomogeneous 
wave  equation  (3.17)  can  be  written  as,  Ez  =  Ez  g  +  Ezp,  the  sum  of  the  general  solution  EZ}9  of 
the  homogeneous  wave  equation  (i.e.,  (3.17)  when  setting  f(z,t )  =  0)  and  the  particular  solution 
Ez^p,  [Evans,  2010],  where 


Ez,P  = 


yfW 


t  rz+(t.-t')  / 


0  J z—(t—t')/ y/JIe 


f(z,  t)dzcLt 


(3.18) 


Eqn.  (3.18)  reveals  that  the  electric  field  computed  in  the  forward  model  depends  on  the  excitation 
source,  which  also  affects  the  reconstruction  procedure  for  tomographic  imaging.  However,  an  exact 
relationship  between  the  excitation  sources  Jz  and  the  quality  of  imaging  e(r)  by  iterative  algorithm 
(i.e.,  convergence  speed  and  reconstruction  error)  is  implicit.  Here  we  rely  on  numerical  means  to 
study  how  source  encoding  schemes  effect  the  quality  of  tomographic  imaging. 


3.4  Numerical  Experiments 

In  this  section,  we  conduct  numerical  experiments  to  test  and  verify  our  results.  The  test  example  is 
shown  in  Fig.  3.1.  The  object  to  be  reconstructed  is  a  two-layer  dish-like  objects.  The  inner  object 
is  of  diameter  of  6  mm  with  a  dielectric  constant  e  =  1.2eo  and  the  outer  object  is  of  diameter  of 
12  mm  with  a  dielectric  constant  e  =  1.5eo-  eo  is  the  vacuum  permittivity  8.85  x  10-12  F/m.  The 
permittivity  value  for  the  surrounding  medium  is  eo-  In  the  simulation,  only  the  permittivity  e(r)  is 
to  be  reconstructed.  The  conductivity  value  is  set  to  zero  cr  =  0,  and  the  magnetic  permeability  is 
chosen  to  be  a  constant  /r  =  Eo  =  47T  x  10 7 H/m.  The  antennas  are  placed  on  the  four  sides  of  the 
imaging  region.  There  arc  a  total  of  90  co-located  transmitters/receivers.  The  computational  area  is 
a  square  region  having  sides  of  length  12  cm  and  is  discretized  to  40  by  40  grids.  In  the  computation, 
we  use  the  MUR  absorbing  boundary  condition  [Mur,  1981].  5%  of  random  disturbance  as  rjm  is 
introduced  to  the  simulated  electric  wavefield  signals.  To  evaluate  the  performance  of  the  iterative 
algorithms  with  different  encoding  schemes,  we  calculate  the  relative  error 

(3  =  ||e(fc)  -  etrue||2/||e(0)  -  etme||2  (3.19) 
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Figure  3.2:  Reconstructed  images  of  dielectric  value  e.  Top:  random  phase  encoding.  Middle: 
waveform  delay  encoding.  Bottom:  uniform  weight  encoding. 


where  and  e('k>  is  the  initial  and  final  value  of  the  reconstructed  image,  respectively.  etme  is  the 
true  value  of  the  dielectric  constant  of  the  target  and  the  medium. 

Simulations  are  conducted  using  Lm  =  7  excitation  sources  at  each  transmission  under  the 
three  different  waveform  encoding  schemes.  Fig.  3.2  depicts  three  reconstructed  images  under 
the  random  phase  encoding  (3.13),  waveform  delay  encoding  (3.14),  and  uniform  encoding  (3.15). 
Fig.  3.3  depicts  the  corresponding  convergence  history  (i.e.,  relative  residual  error  3  vs.  the  CPU 
run  time  in  seconds)  of  the  three  iterative  algorithms.  The  conventional  single  source  excitation 
scheme  is  also  included.  Among  the  three  encoding  schemes,  the  random  phase  encoding  achieves 
fastest  convergence  and  smallest  imaging  error,  compared  with  the  uniform  waveform  encoding 
and  the  waveform  delay  encoding  schemes.  As  expected,  the  conventional  SIMO  imaging  yields 
the  slowest  convergence  and  lowest  image  quality  at  the  450  seconds  compared  with  other  three 
encoding  techniques,  as  shown  in  Fig.  3.3. 

3.5  Conclusions  and  Discussions 

We  developed  waveform  encoding  techniques  for  a  nonlinear  electromagnetic  tomographic  imag¬ 
ing  problem  in  the  time  domain.  Proper  encoding  techniques  such  as  random  phase  encoding  yields 
fast  convergence  of  iterative  reconstruction  while  achieving  high  quality  images  with  reduced  data 
volumes,  thus  demonstrating  the  power  of  signal  processing  techniques  for  improving  significantly 
the  computational  efficiency  for  solving  nonlinear  inverse  problems.  We  note  that  randomized  re¬ 
construction  for  nonlinear  imaging  when  the  reconstruction  step  is  optimized  with  regard  to  each 
random  weights  remains  an  open  research  problem.  The  theory  of  compressive  sensing  for  linear 
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Figure  3.3:  Convergence  plots  of  iterative  algorithms:  SIMO  vs.  three  MIMO  waveform  encoding 
schemes 

inverse  problems  offers  interesting  insights  and  could  lead  to,  combined  with  waveform  encoding 
techniques,  new  fast  and  efficient  EM  full  wave  tomographic  imaging  algorithms  for  large  scale 
imaging  applications. 
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